Nonlinear effective medium theory of disordered spring networks 
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Disordered soft materials, such as fibrous networks in biological contexts exhibit a nonlinear elastic 
response. We study such nonlinear behavior with a minimal model for networks on lattice geome- 
tries with simple Hookian elements with disordered spring constant. By developing a mean-field 
approach to calculate the differential elastic bulk modulus for the macroscopic network response of 
such networks under large isotropic deformations, we provide insight into the origins of the strain 
stiffening and softening behavior of these systems. We find that the nonlinear mechanics depends 
only weakly on the lattice geometry and is governed by the average network connectivity. In partic- 
ular, the nonlinear response is controlled by the isostatic connectivity, which depends strongly on 
the applied strain. Our predictions for the strain dependence of the isostatic point as well as the 
strain-dependent differential bulk modulus agree well with numerical results in both two and three 
dimensions. In addition, by using a mapping between the disordered network and a regular network 
with random forces, we calculate the non-affine fluctuations of the deformation field and compare 
it to the numerical results. Finally, we discuss the limitations and implications of the developed 
theory. 
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PACS numbers: 



I. INTRODUCTION 



Rich elastic behavior is a common feature of many soft 
materials such as foams, granular packings and soft 
glasses P, Q , as well as networks of protein fibers that 
form major structural components of cells and tissue 
One characteristic these varied systems share is 
their particular sensitivity to external stress; in densely 
jammed systems, for instance, the external pressure can 
cause the system to transition between rigid and floppy 
states P, while reconstituted biological filamen- 

tous networks exhibit dramatic strain stiffening under 
shear 

EMU- 

This remarkable nonlinear elastic behav- 
ior of fiber networks has attracted much attention in the 
last decade; in addition to the physiological relevance of 
this nonlinear elastic response for cells and many bio- 
logical tissues (Til [l5| , these systems are also interesting 
from a fundamental perspective, owin g to their unusual 
nonlinear materials properties [FlI4l3l , [Ra - l24| , including 
negative normal stresses (25|. Understanding how their 
intrinsic disordered nature affects the elastic deforma- 
tions is required for a complete theoretical description 
of their nonlinear mechanical behavior. Although struc- 
tural disorder and inhomogeneous deformations clearly 
play a central role in jamming systems [U, @, [26j|, their 
precise role in the nonlinear behavior of fibrous networks 
remains unclear @, QJ, d M, S3, SI • 

Prior work on the nonlinear elasticity of random spring 
networks has focussed on triangular lattices with internal 
stresses (2^ ■ In the limit of small disorder, a perturbation 
theory was applied to describe the nonlinear elastic re- 
sponse of such systems with small dilution. It was shown 
numerically that the transition value of the mean coor- 
dination number, at which the network acquires rigidity, 
shifts with applied strain. Interestingly, the perturbation 



theory also appeared to capture the behavior observed 
numerically even for highly diluted networks, since the 
bulk modulus was found to increase linearly with the 
mean coordination number beyond the rigidity perco- 
lation point. Recently, similar nonlinear behavior was 
analysed for random spring networks in jammed config- 
urations. Consistent with prior work on triangular lat- 
tices [2!|, it was shown that this nonlinear response is 
controlled by the central force isostatic point p}; this 
isostatic point characterizes the average connectivity, z, 
at which the number of central-force constraints balances 
the number of degrees of freedom in the system and is 
given by zq = 2d in d dimensions [30| . For jammed sys- 
tems, it was shown that the nonlinear response close to 
this isostatic point is well described by a mean field scal- 
ing approach 0- By contrast, the systems we consider 
here are not in jammed configurations, but instead fall 
in the class of lattice-based rigidity percolation problems 
[3ll - [33l |. For instance, the linear elastic response of fiber 
networks is also governed by the central-force isostatic 
point for a broad range of network connectivities — even 
with fibers with non-central fiber bending interactions — 
but with non-mean-field behavior (34|. Motivated by 
these results, we investigate the nonlinear behavior of 
fiber networks in the limit of vanishing bending rigidity 
under large deformations. 

Here we investigate the nonlinear elastic response 
of random spring networks under isotropic expan- 
sion or compression over a broad range of network 
connectivities — both above and below the small strain 
isostatic point Zq , as illustrated in Fig. [T] From simula- 
tions we find that disordered subisostatic spring networks 
exhibit significant strain-stiffening. We gain insight into 
the origins of this behavior by developing an effective 
medium theory (EM theory) for the nonlinear responses 
of random spring networks on lattice geometries. The 
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nonlinear behavior of these systems depends only weakly 
on network geometry and appears to be controlled largely 
by the mean network connectivity, z. and the applied 
strain, e. Within the framework of this central-force 
network model, the network's stiffness exhibits a tran- 
sition on the two-dimensional phase diagram in e and z, 
which characterizes the strain dependence of the isostatic 
point, a shown in Fig. O The transition curve, z c (e) — 
representing the transition between a rigid and a floppy 
state — is derived using the EM theory and found to agree 
well with the numerical results. Interestingly, the mean- 
field solution predicts that a superiso static central-force 
network looses rigidity and collapses beyond a thresh- 
old in compressional strain, as was observed for perfect 
two-dimensional triangular network in Monte- Carlo sim- 
ulations at low temperatures jH, H(|. The theoretical 
predictions are verified using numerical calculation. 

The mentioned above results for random spring net- 
works may lend insight in the mechanics of biopolymer 
networks. The nonlinear elasticity of reconstituted net- 
works of intracellular biopolymers such as filamentous 
actin (F-actin) and intermediate filaments has in many 
cases been accounted for by the affine entropic model 
(ill , [l2l , l37l . [Hj]. In this model, network disorder is ig- 
nored by assuming a uniform (affine) deformation and, 
consequently, the nonlinear network response is directly 
determined by the nonlinear entropic force-extension be- 
havior of the individual filaments. By contrast, there 
is increasing evidence that the strain-stiffening behavior 
of networks consisting of stiff thick fibers, such as col- 
lagen and bundled actin networks is governed by collec- 
tive non-affine fiber bending deformations (T3. l28ll39ll40j . 
Despite extensive analytical and numerical investigations 
@, El EI IM Hi H3, gJEI, the principles of such net- 
work deformations and the resulting nonlinear network 
response are still unclear. To investigate the implications 
of our results on random spring networks for biological 
filamentous networks, we include a finite bending rigidity 
for the fibers in our network model. For, sufficiently small 
bending rigidities, these networks exhibit nonlinear elas- 
tic behavior. However, in this case the nonlinear behavior 
is not due to a transition between a floppy and a rigid 
state, but between a soft bending-dominated elastic be- 
havior and a stiffer stretching-dominated behavior [l6| . 
Both with and without bending rigidity, the nonlinear 
response is still governed by the central force isostatic 
point. 

The outline of the paper is as follows. In Sec. [XT] we 
define the model and summarize the applied approach 
and the main results of the paper. In Sec. IHII we present 
the mean-field approach in detail, derive the differential 
bulk modulus of a system, analyze the non-affine fluctu- 
ations and, by using a self-consistency check, we identify 
the range of applicability of the performed approxima- 
tions. In Sec. IIVI we demonstrate the presented general 
method using a particular example of the diluted regu- 
lar networks and compare the analytical predictions to 
the numerical results. We discuss the results and their 




Figure 1: (Color online) A small section of the relaxed 
(left) and expanded (right) diluted triangular lattice. The 
average coordination number in this example is z = 
3. On the right plot the blue/red colors mixture repre- 
sents the high/low elastic energy stored in a bond. See 
http: //www. youtube . com/watch? v=ANSQePygf YU for a stiff- 
ening movie. 
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Figure 2: (Color online) The schematic phase diagram for the 
rigidity of random spring networks under an isotropic strain e. 
The central- force isostatic point, zo, the conductivity thresh- 
old, z COI1 d and the lower rigidity threshold in negative strain 
are indicated in the diagram. 

implications and summarize in Sec. IVI1 

II. THE MODEL 

In this paper we analyze the nonlinear elastic behavior 
of a random central-force network on lattice-geometries 
with varying connectivity. We start out with a model 
in which the bending energy of the fibers or bonds is ig- 
nored. This model will allow us to study the effects of 
finite strain on the central- force isostatic point, z c , and 
the stretching energies of the bonds. To further sim- 
plify our model we only consider isotropic expansional 
and compressional strains (see Fig. [TJ. The calculation 
of the elastic properties under nonlinear shear is compli- 
cated by the broken symmetry and is described elsewhere 
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[47| . The network is constructed on an ordered lattice ge- 
ometry in d > 1 dimensions. We capture the effects of 
disorder by a distribution of the spring constants associ- 
ated to the bonds in the network. In this model the rest 
lengths of all springs are chosen to be identical and equal 
to the lattice spacing, Iq. 

Measuring all lengths in units of £q, the Hamiltonian 
of the system is given by 



2< 



H = -2jHj 



u, 



I) 2 



(1) 



(ij) 



where is the position of vertex i, (ij) denotes the sum- 
mation over neighbouring lattice vertices and /iy is the 
stretching modulus for the bond between vertices i and 
j. The stretching moduli /Uy are drawn independently 
from a known probability density, P (/iy ) . 



A. Quantities of interest 

Here, we investigate the elastic response of the network 
to an applied global expansion/ compression strain, 



L'-L 



(2) 



where V and L are the linear size of the strained and 
unstrained networks, respectively. To quantify the non- 
linear elastic response to the global bulk strain, we define 
the nonlinear differential bulk modulus as 



B(e) 



n d 2 E (e) 
d 2 Qe 2 ' 



(3) 



where E (e) is the average elastic energy per bond, n 
is the density of bonds in the unstrained and undiluted 
lattice and d is dimension of the system. For example, 
for the FCC lattice n = 2v / 6^ 3 , while for the triangular 
lattice n = 2y^3£g 2 . 

This definition of the nonlinear bulk modulus has the 
following advantages: 

a. Other quantities, related to the nonlinear re- 
sponse of the system to a global strain may be deduced 
from B (e) . For instance, the pressure, 



n 



dU 



(4) 



where U = NE is the total elastic energy, N is the to- 
tal number of springs in the undiluted network, V = 
Vq (1 + e) d is the system's volume and Vq is the total 
volume of the unstrained network. This pressure can 
be obtained directly from the nonlinear differential bulk 
modulus using 



n 



(l + e)° 



- f B(e)de. 
Jo 



(5) 



b. In the linear regime, e — > 0, the nonlinear bulk 
modulus converges to 



B(e^ 0) = V 



d 2 u 

dV 2 ' 



(6) 



c. If the material is composed of Hookian bonds 
and its deformation is affine, B (e) is constant and equal 
to n/d 2 times the average spring constant of the network. 
Thus, by plotting =-B (e) one can easily compare the 
actual elastic properties to the affine predictions. 



B. Numerical results 

To study the nonlinear elastic response of random 
spring networks, we choose a specific realization of the 
spring constant probability density for networks on lat- 
tice geometries (see Sec. [IV] for a detailed discussion). In 
particular, we investigate bond diluted network of springs 
with a modulus [i with the following probability density 



P (jHj) = VS {puj -v) + (l-V)d Oy) 



(7) 



Thus, either a bond with a stretch modulus [i is present 
with a probability V ', or absent with a probability 1 — V . 
By varying V we can tune the connectivity, i.e. the mean 
number of springs, z, which are attached to a crosslink of 
the network. Here we study nonlinear elasticity of such 
bond-diluted networks on a triangular lattice in d = 2 
and face centered cubic (FCC) lattice in d = 3. 

The mechanical response of these networks is sensitive 
to the applied strain. We quantify this network response 
with a differential bulk modulus, B (z, e) as shown in 
Fig. [31 The qualitative behavior of the nonlinear bulk 
modulus is similar to the behavior of the bulk moduli in 
the linear regime, but with an isostatic point that shifts 
continuously to lower coordination numbers, as demon- 
strated in Fig. [3] The nonlinear response of the di- 
luted central-force networks can be characterized by a 
two-dimensional (z, e) phase diagram of the differential 
bulk modulus, as shown schematically in Fig. [5J The 
transition curve, z c (z), separates the floppy and the rigid 
regions. From this it can be understood that a subiso- 
static diluted regular network with central-force interac- 
tions exhibits a strain-stiffening behavior — from a floppy 
to a rigid structure — as a function of applied strain. 



C. Mean-field approach 

We gain insight in the nonlinear response of random 
spring networks (see Eq. {1])), by developing an effec- 
tive medium approach for the high strain regime. This 
EM theory aims to provide a complete quantitative de- 
scription of the nonlinear elastic response of a network 
under an external expansional/compressional strain e. 
Our approach is a nonlinear extension of the EM ap- 
proaches used to successfully describe the linear elastic 
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Figure 3: (Color online) The nonlinear differential bulk mod- 
ulus as a function of the coordination number for triangular 
(open symbols) and FCC (filled symbols) lattices for different 
strain values. The solid lines are the results of the EM theory 
(Eq. d33) or, in the explicit form, Eq. (JBTJ). 



response of diluted lattice based networks [48|, and goes 
beyond perturbative approaches for networks with small 
dilution [2§|. The effective medium theory for the linear 
elastic response of the disordered spring networks under 
small deformation was shown to predict the location of 
the critical coordination number and the elastic response 
far from it [48l450l | . In other systems with non-central 
force interactions, such as fiber bending models, the EM 
theory succeeded to capture the qualitative elastic be- 
havior of the network [34j. 

The nonlinear EM approach developed here is based 
on a scheme to construct a mapping from the disordered 
system onto a perfect lattice system with uniform bond 
stiffness with the same network topology and strain e, as 
illustrated in Fig. 2J This mapping is realized by an effec- 
tive uniform central force interaction, /iy — > Jl. The effec- 
tive parameter, Jl, is determined using a self-consistency 
requirement: replacing a random bond in the uniform 
EM under strain with a bond drawn from the original 
probability density, P (^ij), results in a local fluctuation 
in the deformation field, which vanishes when averaged 
over the distribution P (nij). In addition, we assume that 
the fluctuations of the deformations are small compared 
to the distance between crosslinks. Importantly, the ef- 
fective parameter Jl depends on strain, which gives rise 
to a nonlinear elastic response. 



III. EFFECTIVE MEDIUM THEORY 

We apply the EM theory method to a network sub- 
jected to a global expansion with a macroscopic isotropic 
strain e. The position of a crosslink i is given by 



strained configuration and is the displacement due to 
the applied strain. The affine displacement is defined as 
uf ff — u^ ff = (1 + e) Tij, where iy,- is the vector from R" 
to R° in the undeformed reference state. Here we allow 
for non-affine displacements 



V,; = U, - U 



ail' 



(8) 



given by the deviation of the displacement of network 
node i from its affine value. However, we assume that 
the resulting non-affine relative displacements of neigh- 
bouring nodes i and j are much smaller than the corre- 
sponding affine displacement: 



= v,- — v, 



< |uf ff - uf\ 



t e. 



(9) 



Thus, we can expand the Hamiltonian around the affine 
strain configuration (small i>ij). Up to second order in 



Vij we arrive at [2£ 



2 e +ev *- 



1 (vy • r ;i ) +1 

2 l + e 



(10) 

The expansion of the whole network corresponds to the 
global constraint 



5> y = o. 



(ii) 



R, 



RV 



Ui, where R" is the position in the un- 



interestingly, the Hamiltonian for the network under 
a finite strain bears a resemblance with the Born model, 
which includes both isotropic and anisotropic pairwise 
interactions 1 5 ill — as is the case here. However, as was 
noticed in [29j, the model presented here is formally 
distinct from the Born model, since here is the dis- 
placement from the affinely deformed state and not from 
the undeformed configuration as is the case in the Born 
model. Thus, from our model (Eq. <jT0j) ) we see that 
a finite strain introduces additional interactions that pe- 
nalize non-affine deformations with a coupling parameter 
that is directly proportional to the strain e. 
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To investigate the nonlinear elastic behavior of the 
model in Eq. (| 10[) , we set up an effective medium theory. 
In the EM approach, we mimic the disordered system by 
the regular one with an effective parameter, i.e. faj — > Jl. 
To deform the EM network similarly to the original sys- 
tem with the same global strain e, one can use a Lagrange 
multiplier to ensure that the global constraint (fTTj) is sat- 
isfied. In other words, the EM network may be globally 
expanded by applying the force that assures mechanical 
equilibrium for the affine, = 0, configuration. Thus, 
the EM system has the Hamiltonian, given by 



Hem 



(ij) 



1 (Vij 

2 



(12) 



where Ay = — juery. To calculate the effective param- 
eter Jl we demand self-consistency of the EM j48|. The 
self-consistency requirement in this context can be for- 
mulated as follows: the non-affine displacement induced 
by the replacement of a single bond in the EM vanishes 
on average, 



(v,. 



= 0. 



(13) 



Here, the average is taken over the distribution of the 
nm bond in the original disordered system, i.e. accord- 
ing to the probability density P (u nm )- To calculate the 
displacement v„ m after the replacement we solve the per- 
turbed EM Hamiltonian that is given by 



Hem +^ 0" 

nm 



8 1 



,) 2 +ev, 



2 

lit n 



(14) 



In the configuration that minimizes the energy, the dis- 
placement of the nm bond is given by 



UEM + fJ-nm — A 4 ' 



(15) 



where uem is the displacement of the nm bond in the 
unperturbed EM network due to a unit force acting along 
the nm bond. As detailed in Appendix fAJ it is given by 



Hem 



a(eY 



(16) 



where a (e) is given in Eq. (|A6[) and may be approxi- 
mated for a highly coordinated lattice by 



o(e) 



2d(l 



'-3 



1 



d-l 



2+d 



+ £ 



(17) 



to the following equation for the effective parameter 1 



1 - Jl (e) I Hi 



^ + 1 - fi (e) I m 



-P {Hij)duij = 0. 



(18) 



Importantly, in contrast to the linear EM, here the 
effective parameter Jl (e) cannot be interpreted as the ef- 
fective spring-constant of the bonds in the EM. However, 
using the expression for Jl (e) one can determine the elas- 
tic properties of the original disordered system as follows. 
Since the equilibrium configuration of the regular, EM 
network is given by the affine expansion, vy = 0, its 
energy (per bond) is given by 

Eem (e) = Hem ^ = 0) = ^ (e) e 2 . (19) 

The last expression may be interpreted as an approxi- 
mation for the original system's energy up to correction 
terms. These terms appear since the energy is defined as 



E EM (e) = - f f B EM (e")de'' 
n Jo Jo 



or 



E 



EM 



(e) = 



n 



EM 



(e')e'. 



(20) 



(21) 



Thus, for z < zq it includes the integration in the floppy 
phase, where the EM theory breaks down and predicts 
non-physical elastic properties. To calculate the correc- 
tion terms we calculate first the floppy and rigid phases 
separation curve and then subtract from the expression 
(UHJ) the integration in the floppy phase. For a given 
mean coordination number the transition strain value, 
e c (z), may be found as follows. If the transition is first 
order, one requires that the energy and its first derivative 
vanish at the transition point. By contrast, if the tran- 
sition is second order, one requires that the energy and 
its first and second derivatives vanish at the transition 
point. These two possible assumptions about the transi- 
tion order result in different transition curves (z Cl e c ). We 
define them as (z Cl , e Cl ) for the first order transition case 
and (z C2 ,e C2 ) for the second order transition case. Since 
the order of the transition cannot be deduced from the 
EM theory, described here, we will analyze both options. 
In Sec. IIVI we calculate both transition curves for the 
particular example of the diluted regular networks and 
compare them to the numerical results. 

The nonlinear bulk modulus, defined in Eq. ([3]), does 
not depend on the transition order and may be approxi- 
mated by the EM approach using Eqs. (flU)) and (|2TJ)) . It 
is given by 



R (A - U 9 



M(e) y 



(22) 



Given Eqs. (|15I16|) . the self-consistency Eq. (fT3"|) leads 1 In the linear regime, Eq. ([T8j reduces to Eq. (9) in Ref. (48 
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where Jl (e) is given in Eq. (|A6[) and is approximated 
in Eq. (|17p . In the limit of small strain, given by 
e —> 0, the mean-field rigidity threshold of the network, 
z c (e — > 0) = 2d, also does not depend on the transition 
order. This corresponds to the mean-field isostatic point 
(30l . HH]. In the limit of large strain, e — > 00, the mean- 
field rigidity threshold of the network does not depend 
on the transition order and is z c (e — > 00) = 2. 

The proper energy and the pressure in the system can 
be obtained from Bem (e) by integration from e c 



d 2 f e r' 
Eem (e) — — / / B E m (e") de"de 



(23) 



The obtained regular lattice (EM) with random forces is 
equivalent to model B in Ref. (52| (in contrast to the orig- 
inal system which is defined as model A in Ref. (52|) and 
may thus be treated similarly. In particular, the non- 
affinity correlation function of the strain field may be 
evaluated within this model. The described mapping has 
the same level of approximation as the EM approxima- 
tion. However, this mapping completes the EM theory 
in the sense that it allows to calculate all the correla- 
tion functions. Below we provide two examples of the 
usefulness of this mapping by calculating two important 
one-point correlation functions: the average non-affine 
displacements and its nonlinear, differential analog. 



and 



n 



EM 



(1 



B 



EM 



(e')de', (24) 



where e c is defined as e Cl or e C2 for the first and second 
order transition assumptions, respectively. 

The approach presented in this section allows one to 
calculate the elastic parameters of a system with a given 
topology and elastic constant distribution in the nonlin- 
ear elastic regime. The nonlinear differential bulk mod- 
ulus is given by Eq. (|22[) while Eq. (fT5|) determines the 
effective parameter Jl(e). Eq. Q18p may be solved nu- 
merically for any realization of the spring constant prob- 
ability density, P(fiij). In Sec. IIVI we demonstrate the 
presented method using the particular example of diluted 
regular networks when Eq. (jI8j) can be easily solved an- 
alytically. 



A. Mapping fluctuations to random forces 

Within the EM approach, desribed above, a deviation 
of the spring constant of a given bond from the EM spring 
constant is described as an additional force dipole that 
act on this bond. Therefore, the disorder of the spring 
constant may be mapped to the disorder of the force 
dipoles which act on the regular EM in the EM theory 
approximation. More specifically, the replacement of the 
EM spring between nodes i and j by the spring with the 
elastic constant fiij is equivalent to the force 



fi - fiij 



/' 



a(e) V + I 1 *! 



a(e) 



(25) 



acting along the bond. Due to the self-consistence re- 
quirement (|13p the average force is zero, 



<tf>=0 



(26) 



and, since we assumed that there is no correlation be- 
tween spring constants on distinct bonds, the associated 
forces are also uncorrelated: 



(27) 



B. Correlation functions 

It is known that the rigidity percolation transition in 
the small-strain limit is accompanied by highly non-affine 

network deformations @, m, m. As we show in this 

work, the non-affine deformation field is responsible for 
the strain stiffening behavior of the disordered spring net- 
works. The mapping, described in Sec. MI Al allows us 
to calculate any correlation function of the displacement 
field. 



1. N on- affinity parameter 

Several methods have been proposed to quantify the 
deviation from a uniform (affine) strain field pjl |42|, [52|, 
[53| . A useful parameter for the non-affinity characteriza- 
tion is the average deviation from the affine configuration, 



(28) 



where (••■}„ denotes the average over all vertices. Be- 
low we calculate r using the EM approach including the 
mapping described in Sec. IHI Al 

The Fourier transform of the force (f2"5")) is given by 



P 3 ( k ) = e M _t_ (gik-R, _ e ik-R^ (2Q) 



Sft - A* + Mij 



a(e) 



Thus, the Fourier transform of the displacement field 
from the affine configuration due to this force is given 
by 

v(k) = -^D-^f^k) 

m 



e (1 + e)^ ^T-M+M»j 



f-f £(r®r + 6l)(I - e^ T ) 

W> r 



■ * 1 
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or in real space 

_£(!+£) 
V " Na(e) 



E 



(r <g> r + el) (1 - e' lk ' r ) 



(31) 

Since the forces are independent and identically dis- 
tributed random variables, the variance of the displace- 
ment from the affine configuration of every network 
crosslink is given by 



(v 2 ) =- 



e(l + e) 




(J, - fJLij 



sfa - v- + , 



-T 



1 - e l 



X) (r ® r + el) (1 - e ik ' r ) 



(32) 



The same approximation of highly coordinated lattice 
that was used to derive Eq. (|A7[I now gives 



2Za 2 (e) 



2 + d 




(33) 



The sum over k may be estimated for any dimension 

k 



where 




(35) 



d > 2 



and 



A* 



IV i 

k 



0(1) (36) 



is a dimensionless parameter which depends weakly on 
the lattice geometry (for instance, — 0.36 for the tri- 
angular lattice). As defined above, £o is the rest length 
of a bond (throughout this paper this quantity was set to 
unity, but is shown here explicitly to emphasize the cut- 
off of the integral in Fourier space and the units of the 
non-affine parameter) and L = IqN 1 ^ is the size of the 
unstrained network. In sum, the non-affinity parameter 
is given by 



r(e)~ A 



2Za 2 (e) 



3 

2+d 



d-1 

2+d 



a(e) 



fl + fj,ij 



3 

2+d 



1 

2+d 



fd(N) 



(37) 



Differential non-affinity parameter 



Eq. (|3Z 



In the nonlinear regime the more interesting quantity 
is the differential non-affinity fluctuation defined as 



8T(e) = 



v„ (e + Ae) - ±^±^v„ (e) 




r(e). 
(38) 

The first term in the last expression is calculated in Ap- 
pendix [C] while the last two may be easily deduced from 



C. Ginzburg criterion 

Although the mean-field approach is not a controlled 
approximation and has no small parameter, one can 
check for self-consistency of the assumption of small fluc- 
tuations In our case we assume small relative de- 
viations of the EM from the affine strain field (see Eq. 
(HOD): 



(v 2 ) 



(rim) 



< i, 



(39) 



where (•••)/ nm \ is the average over all connected nodes of 
the EM network. Therefore, it is instructive to analyze 
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the behavior of the two point, nearest neighbour corre- 
lation function. Using the same mapping to the random 
forces model as above one gets 



(v 2 ) 

\ nm I 



NN 



(nm) 



r B d 

A d f d d 2 



(40) 



where 



B d 



1 / e o h d-li.2 



/L 



k 2 dk 



Si 



1/e ° k d -Hk 



(41) 



k Jl/L 

The nonlinear, differential version of Tjvjv, defined as 

l+e+Ae, 



ST 



NN 



lim 

At->0 



v nm (e- + Ae) - 



1+6 



Ac 



is given by 



ST 



NN 



ST 

A d f c 



B d . 



(nm) 

(42) 



(43) 



In Sec. HVDI we analyze the non-affinity parameters for 
the particular example of the diluted regular networks 
and compare the analytical results with the numerical 
simulations. 



where (■■•), denotes an average over all network vertices 
and Z is the coordination number given that all exist- 
ing springs have a non-zero spring constant. Below the 
so-called isostatic threshold of the average coordination 
number, zq = z c (e — > 0), the network is floppy (30l . \5l\ . 

For the unstressed reference state and zero strain limit, 
Maxwell introduced a mean-field counting argument for 
this threshold coordination number at which the number 
of degrees of freedom and the number of constraints due 
to the central force interactions are equal. This yields an 
EM approximation for the isostatic coordination number 



2d. 



(45) 



It was conjectured (see Refs. in [48J) that the bulk mod- 
ulus of the diluted network in the zero strain limit can 
be expressed in term of zq as 



B (e -> 0) = ii 



d 2 Z 



(46) 



Eqs. (|45l46p were derived [48| using the EM theory and 
were shown to predict well the location of the isostatic 
point and the elasticity of a diluted network far from its 
isostatic point. 



B. Infinite strain limit 



IV. DILUTED REGULAR NETWORKS 

A significant understanding of different physical phe- 
nomena in disordered systems, including percolation 
55 , l56l | and the elastic behavior of amorphous materials 
56| was achieved by modelling the topological disorder 
by a random dilution of a regular structure. Motivated 
by this we demonstrate the mean-field solution presented 
above using the particular example of bond-diluted regu- 
lar networks. The probability density for the spring con- 
stants for such a network is given in Eq. ([7]). Networks 
of this kind are referred to as diluted spring networks or 
the central-force elastic percolation model. The linear 
elastic response of this model has been extensively stud- 
ied {32I . Here we show how these results generalize 
for large strain values. Before presenting the full mean- 
field solution for these networks, below we briefly sum 
up the known relevant results in the small strain regime 
and discuss the infinite strain limit expectations from the 
nonlinear EM theory. 



Before we turn to the full problem with arbitrary strain 
values, it is instructive to discuss our expectations in the 
infinite strain limit. In this limit, the rigidity threshold 
(isostatic point) can be expected to approach the conduc- 
tivity threshold, denoted here by z cond . By analogy with 
the behavior at small strains Eq. (|46j) , we anticipate (and 
derive this result below) that in the large strain limit the 
nonlinear bulk modulus is equal to 



B (e — > 00) = fx- 



(e 



ft Z Z CO nd 



d 2 Z - z c (e -> 00) (PZ- Zcond 

_ (47) 
The mean-field calculation [58| of the conductivity per- 
colation and our calculation below in the infinite strain 
limit both suggest z c (e ^ 00) = 2. We expect a devia- 
tion of the EM theory from the numerical calculation in 
the infinite strain limit close to the conductivity percola- 
tion point due to the failure of the mean-field approach 
to predict the precise value of the conductivity threshold. 



Full solution 



A. Zero strain limit 

The average coordination number for bond diluted net- 
works is defined as 



£(i-<Wo)) =*p, 



(44) 



Using Eq. (7} , the solution for the self consistent equa- 
tion (fT5|) is given by 



H(e) = fi- 



a(e)Z 



a(e)Z' 



(48) 



where a (e) is given in Eq. (|A6[) and is approximated 
in Eq. p7p . The nonlinear differential bulk modulus 
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Figure 5: The nonlinear bulk modulus for the diluted FCC lattice as a function of the applied strain for different values of the 
average coordination number shown as labels next to the curves. The numerical data is depicted as symbols and the nonlinear 
EM theory predictions are shown as solid lines. Big crosses represent the location of the first order transition as predicted by 
Eq. jB2j 



approach using Eqs. 



and 



0.05 



and is given by 




«j| g 0.02 



_ Y 

~ 1-* \04 

/ ' 



10 



10 1 



Figure 6: The nonlinear bulk modulus for the diluted FCC 
lattice as a function of the applied strain for different values 
of the average coordination number shown as labels next to 
the curves. The numerical data is depicted as symbols and 
the nonlinear EM theory predictions are shown as solid lines. 
For z < 2 the EM theory predicts zero bulk modulus for any 
strain. 



B E m (z, e) 



' n a 2 



-a(e)Z 



Z-a(e)Z 2 



e > e c 
e < e c 



(49) 



of the original system may be approximated by the EM 



where e c is defined as e Cl or e C2 for the first and second 
order transition assumptions, respectively. In Appendix 
[5] we present the explicit result for the nonlinear differ- 
ential bulk modulus and the transition curves (for both 
assumptions for the transition order), based on Eq. ([35]). 

A comparison between this analytic prediction and the 
numerical results is shown in Figs. I5I6I7I Below the con- 
ductivity percolation threshold, z < z con d, a network, 
does not resist deformation for any strain and B (e) = 
due to the absence of an infinite connected cluster. By 
contrast, when the coordination number is in the range 
Zcond < z < zq, a network only develops a non-zero differ- 
ential bulk modulus for positive strains above a threshold 
e c (z). For superisostatic coordination numbers, z > zq, 
the differential bulk modulus is larger than zero in the 
small strain limit; B increases with e until it reaches a 
plateau of the large strain limit (see right panels in Figs. 

EE). 

In contrast to the positive strains, for a negative values 
of the strain, the modulus B (2, e) of superisostatic net- 
works decreases with |e| until it vanishes below a thresh- 
old, predicted by the nonlinear EM theory in Eq. (|B3[) 
(see left plots in Figs. I5I7I8[) . This collapse was observed 
for perfect two-dimensional triangular network in Monte- 
Carlo simulations at low temperatures (35l. |36|. Here we 
show that reduction of the mean coordination number 
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Figure 7: The nonlinear bulk modulus for the diluted triangular lattice as a function of the applied strain for different values 
of the average coordination number shown as labels next to the curves. The numerical data is depicted as symbols and the 
nonlinear EM theory predictions are shown as solid lines. Big crosses represent the location of the first order transition as 
predicted by Eq. (jBlj) 



shifts this collapse towards smaller values of |e|. 

The agreement with the numerical data is good far 
from the transition point. The bulk modulus in the infi- 
nite strain limit is given by 



B 



EM 



(e — > oo) = /i 



d? Z-2' 



(50) 



such that the transition average coordination number ap- 
proaches the conductivity threshold, z c (e — > oo) = 2. 
However, the mean-field prediction for the conductiv- 
ity threshold deviates from the numerical result [HD, [5j| . 
This may account for the discrepancy between the non- 
linear EM theory prediction and the simulation results 
close to the conductivity percolation threshold in the 
large strain regime (see Fig. [5] and large strain values for 
d = 3 in Fig. [5]). In fact, we find that for the FCC lattice 
the rigidity percolation in the large strain limit occurs at 



z c (e — > oo) = z c 



1.5 ± 0.3. This result is consistent 



with both the empirical law for the conductivity thresh- 
old z con d ~ |59fl and with the numerical result of the 
FCC lattice conductivity threshold z con d ~ 1.442 [60l |. 

The results discussed above are summarized in a phase 
diagram shown in Fig. [5J The curves indicate the transi- 
tion connectivity number between rigid and floppy phases 
as a function of applied strain. The explicit formulas may 
be found in Appendix [B] for both assumptions about the 
transition order (see Eqs. (|T32|) and (|B3[) V The strain 
dependence of the isostatic point we find numerically is 
reasonably well described by the nonlinear EM theory, 
as shown in Fig. [5] For the negative strain values at 
the transition point only the differential bulk modulus 




Figure 8: The phase diagram for triangular (upper panels) 
and FCC (bottom panels) lattices. The numerical data for the 
transition points is depicted as symbols and the nonlinear EM 
theory predictions are shown as solid lines for the first order 
transition given by Eq. (|B2[) and dashed lines for the second 
order transition given by Eq. (|B3[) . The curves separate the 
floppy (below) from the rigid (above) phases. 



vanishes but not the stress and the elastic energy; this 
unambiguously corresponds to a second order transition. 
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D. Non-affine fluctuations 

Here we demonstrate the method presented in Sees. 
MI Al and MI Bl and analyze the non-affine fluctuations 
for the particular case of the diluted regular lattices. Us- 
ing Eqs. (|57|) . ([55)1 and the expression for the effective 
parameter, /1(e), Eq. (|48p . one obtains the expressions 
for the non-affinity parameters T and ST. A comparison 
between the analytical formula and the numerical calcu- 
lation is shown in Fig. [51 For superisostatic networks 
the numerical results agree well with the nonlinear EM 
theory predictions. However, on the transition curve, 
the non-affinity parameter appears to diverge; this di- 
vergence is not captured by the nonlinear EM theory. 
Further insight in this discrepancy between the EM the- 
ory and numerical results over a range of parameters (in- 
cluding the divergences) is gained by analysing two-point 
correlation functions. 

As discussed in Sec. MI C[ the mean-field assumption 
of small fluctuations can be determined self-consistently 
using the two point, nearest neighbour correlation func- 
tion defined in Eq. (|40|) . Based on Eq. (|43|) one may 
calculate the non-affine fluctuations, Tnn, and its dif- 
ferential analog, 5Tnn, for the expanded diluted regular 
networks. The comparison between the theoretical calcu- 
lation of r^vAf and STnn, including the numerical results 
is shown in Fig. [TUJ The agreement between the theoret- 
ical prediction and the numerical data is good when the 
Ginzburg criterion is satisfied, i.e. Tnn *C 1. Clearly, 
close to the transition curve, where the non-affinity pa- 
rameters diverge, one can expect the EM theory to fail 
since the Ginzburg criterion is strongly violated. Note, 
however, that the theoretical prediction for the bulk elas- 
tic properties appears to be reasonable even when the 
Ginzburg criterion is not satisfied (see Figs. 171 and [POl for 
small values of z). 

E. Additional weak interactions: fiber bending 

Many biopolymer networks, including collagen and fib- 
rin networks, have a branched structure with a connec- 
tivity close to three on average [4(|. The rigidity of such 
networks with connectivities below Maxwell's central- 
force isostatic point can be accounted for by the exis- 
tence of additional non-central force interactions such as 
those arising from fiber bending. To analyze the effects 
of the finite fiber bending stiffness on network elasticity 
we generalize the model presented in [34[ to the nonlin- 
ear regime. The resulting Hamiltonian is composed of 
two terms representing the stretching and the bending 
energies: 

H = -^gij (K - u j \ - 1 ' '•'X/ y '--' / - 4 ~ cos %fe) ' 

(ij) (ijk) 

(51) 



where the summation in the bending term extends over 
consecutive bonds along the same fiber and Aflyfc is the 
angle between the ij and the jk bonds. Here, g- t j = 1 for 
uncut bonds (with fj,ij = fi) and gij = for bonds that 
have been cut (with fj,ij = 0). Thus, in this model the 
network cross-links are freely hinging. Various EM theo- 
ries for this model were developed for the linear elastic- 
ity [34l , lolT | . However, the generalization to the nonlinear 
regime seems to be technically challenging. The nonlin- 
ear EM theory described above is used to calculate the 
nonlinear differential bulk modulus in the k = case. 
To analyze the importance of the additional interactions 
we compare our purely central-force, n = 0, analytical 
formula Eqs. (|49IBip for the nonlinear bulk modulus to 
the numerical results for different values of k (see Fig. 
Hip . For small enough values of k the nonlinear bulk 
modulus does not depend on k above the central-force 
isostatic point. By contrast, below the transition strain, 
B (e) approaches a plateau proportional to n, which is 
not captured by the central-force nonlinear EM theory. 
However, the strain at which the network exhibits strain 
stiffening appears to be well approximated by the non- 
linear EM theory prediction. 



V. BRANCHED NETWORKS AND RANDOM 
BOND MODEL 

The nonlinear EM theory predicts an expression for the 
bulk modulus (Eq. (|49[) . and more explicitly in Eq. (JbTJ) ) 
for a regular, diluted elastic network that depends on the 
geometry of the undiluted network only via its coordina- 
tion number, Z. This number sets the maximal possible 
coordination number of the lattice. In some cases, such 
as collagen-I that exhibits a branched structure |4Cj[ , the 
maximal possible coordination number seems to be very 
high such that the probability distribution of z is expo- 
nential. For such a networks it is instructive to consider a 
high Z limit of the expression for the bulk modulus, Eqs. 
(US]) or (|B1|) . More specifically, an effectively off-lattice 
network with an exponentially distributed coordination 
number with a given mean 

z = VZ (52) 

can be constructed by almost total dilution, V — > 0, of 
a highly coordinated undiluted regular lattice, Z — > 00. 
A particularly simple example of this type of network is 
the Random Bond Model where randomly located points 
are randomly connected [62|,[63| (see Fig. [T2"|) . Using this 
limiting procedure the expressions for the bulk modulus 
reduces to 
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Figure 9: Non-affinity one-point correlation functions F (red filled circles) and <5F (blue empty squares) as a function of the 
applied strain, e for different values of the mean coordination number, z (see upper right corner of every plot) on the diluted 
triangular lattice of size 140 x 140. The lines represent the theoretical prediction for Y (red solid line) and ST (blue dashed 
line) . 




Figure 10: Non-affinity two-point correlation functions Fjvat (red filled circles) and STnn (blue empty squares) as a function 
of the applied strain, e for different values of the mean coordination number, z (see upper right corner of every plot) on the 
diluted triangular lattice. The lines represent the theoretical prediction for Fjvjv (red solid line) and SFnn (blue dashed line). 
The thick gray line indicate the value of 1 to compare with Fjvat for the Ginzburg criterion. 



13 




Figure 11: The nonlinear bulk modulus for the triangular lattice as a function of the applied strain for different values of the 
bending modulus (see legend above the plots). The average coordination numbers are 3, 4.02 and 4.8 (see upper left corner on 
each plot). The theoretical curves are given by Eq. (|49[) and, in the explicit form, by Eq. (JET}. Insets show the same plots on 
the semi-log scale. Big crosses represent the location of the first order transition as predicted by Eq. (IB2|) 



rid ( d — 1 

B EM (z,e) = ^lz-2- — 



2(d + l) 



[l + (2 + d)e] d (l+4±2 e ) 



(53) 



where rid is the density of the bonds of the diluted net- 
work. The transition curve is given by 



z c (e) = 2 



1 



2(d + i) 



G 



d + 2 \[l + (2 + d)ef + 



(54) 

These results may be directly applied to the Random 
Bond Model. A comparison between the numerical cal- 
culation of the Random Bond Model bulk modulus and 
Eq. (|53| for d — 2 is shown in Fig. [13] Branched net- 
works with exponential distribution of the local connec- 
tivity are also expected to behave according to Eq. (1ST 



VI. SUMMARY AND DISCUSSION 

Motivated by the rich nonlinear elastic behavior of dis- 
ordered materials we studied random spring networks 
under isotropic strains. We provided a quantitative an- 
alytical theory for the strain-stiffening phenomena that 
is driven by non-uniform (non-afhne) deformation fields 
originating in the network disorder. We considered disor- 
dered networks on lattice topologies of Hookean springs. 
The disorder is introduced with a non-uniform distribu- 
tion of the spring constants for the bonds on a regular 
lattices. The central parameter that characterizes such 
networks is the mean coordination number, z; the thresh- 
old, z = z c , separates a floppy from a rigid phases. How- 



ever, when some fraction of the bonds are under stress, 
the isostatic coordination number can shift continuously 
to lower values [57j • This can be realized, for example, by 
applying a large deformation to the network Q or by in- 
troducing local contractile forces (4(| . It was shown that 
rigidity can be induced by additional stresses or strains 
in networks with connectivities below the zero-strain iso- 
static point, zq = z c (e — > 0) = 2d in d dimensions. As 
a result, a significant strain stiffening is induced as the 
network transitions from the floppy to the rigid phase. 

Here we have developed a nonlinear effective medium 
approach for regular central-force networks with disor- 
dered spring constants to provide insight in such be- 
havior. In this model, we expand the Hamiltonian 
around the affine deformation state for an isotropically 
expanded network. Thus, this theory explicitly accounts 
for non-affine deformations that are small compared to 
the affinely strained unit cell. The main result of the EM 
theory approach is the nonlinear differential bulk mod- 
ulus given by Eq. (|49[) . where the effective parameter 
Jl may be found for a given spring constant probabil- 
ity density of the original network, P(fiij), using Eq. 
(HHJ). We demonstrated that this theory quantitatively 
captures the nonlinear elastic properties of a bond diluted 
network for arbitrary strains far from a transition in two 
and three dimensions. In particular, this theory predicts 
a continuous transition curve for the strain dependent 
isostatic point, varying from z$ = z c (e — ¥ 0) ~ 2d at zero 
strain to the conductivity threshold z con d — 2 in the infi- 
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Figure 12: Illustration of the Random Bond Model network 
in two dimensions. Here the density of the bonds is rid = 1 /3 
and the mean coordination number is z = 2. 




Figure 13: The differential nonlinear bulk modulus as a func- 
tion of the applied strain of the two-dimensional Random 
Bond Model network for different values of the mean coor- 
dination number, z. The dots represent the numerical result 
while the solid lines are based on Eq. (|53|l . 



nite strain limit. The transition at the strain-dependent 
rigidity point is accompanied with divergent strain fluctu- 
ations, reminiscent of critical behavior. We showed how 
the nonlinear EM theory can be used to calculate the cor- 
relation functions associated to such strain fluctuations. 
The two point non-affinity parameter, quantifying the 
relative non-affine deformations of neighboring points, 
can be used to inspect internal consistency (Ginzburg 
criterion) of the nonlinear EM theory approach, which 
breaks down in the vicinity of the (strain dependent) 
rigidity percolation point. 

Application of the EM theory developed here presup- 
poses that the order of the transition is known. From our 
numerical results one cannot rule out the possibility of 
either a first or a second order transition. This remains 
a subject of further study (47| . 

We found that a superisostatic disordered network 
with z > 2d may loose rigidity under positive pressure 
(negative strain values) in two- and three-dimensional 
networks. Similar elastic collapse was found and analysed 
for the perfect triangular lattice (35l . |36| . Here we showed 
that the location of such a collapse depends mostly on 
the network topology via the mean coordination num- 
ber. The mean-field approach developed here is found 
to predict reasonably well the location of the collapse 
and the elastic properties of the network for the negative 
strain values. 

We also investigated the effects of additional weak non- 
central force interactions numerically in the form of fiber 
bending in bond diluted networks. The resulting fiber 
network exhibits a strain stiffening transition from a soft, 
bending dominated regime to a stretching dominated 
regime. Importantly, however, this transition still oc- 
curs at the transition strain predicted by the central-force 
nonlinear EM theory, which quantitatively captures the 
nonlinear elasticity beyond the transition strain. These 
results may lend insight into the nonlinear elasticity of 
biological fiber networks. 

The EM theory expressions for the elasticity behavior 
of random spring networks depends on the network geom- 
etry only via the coordination number of the undiluted 
network, Z. One may interpret Z also as the maximal 
coordination number of the diluted network. This may 
lead to the temptation to use the results obtained here 
for other than diluted regular network systems, including 
networks with geometrical disorder. However, there is at 
least one example where a network based on the jammed 
configuration geometry has qualitatively different elastic 
behavior than the diluted regular network with the same 
mean coordination number [64|. In the present work we 
defined the strain of the diluted network relative to the 
zero energy state of the undiluted network. Therefore, 
it is unclear how our results may be extrapolated for ge- 
ometrically disordered networks. Nevertheless, we have 
shown that our results may be applied to describe the 
elastic response of the geometrically disordered Random 
Bond Model. 

In this work we focused on the differential bulk mod- 
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Appendix A: The calculation of (j,em 

In this Appendix we calculate ft em — the displacement of the nm bond in the unperturbed EM network due to a 
unit force r nm acting on the nm bond. 

The dynamical matrix of the unperturbed EM Hamiltonian (|12[) is given by 



Di 



-ifi (r« ® r y + el) i^j 



1 r y + el) i = j 



where I is the unit tensor and (g> is the external product. The Fourier transform of D is given by 

D(k) = ^/V' kr = 



(Al) 



ft \ - 

1 + e^ 



r (X> r + e 



I) (1 



(A2) 



where r runs over all unit bond vectors. The unit force acting on the nm bond is given by 
so that its Fourier transform is 

f (k) = ^J>' k I( = r nm (1 - e ik ' r »-) . 



(A3) 



(A4) 



Thus the Fourier transform of the displacement field is given by 

v(k) = -i?- 1 (k).f(k) 
The displacement of the nm bond due to the unit force is 

1 2d(l + e) 



ftEM 



r nm ■ £V (k) (e-*-'- - 1) = -£> nm ■ f (k) D~ l (k) (e 



-?kr„ 



ft 



1 - 



E(i 



X;(r<8ir + el)(l-e*- r ) 



1) 



o(e) 



(A5) 



(A6) 



For a highly coordinated lattice the sum over r may be well approximated by the integral over the sphere that includes 
all the neighbouring crosslinks and, since the sum over k is dominated by the small k • r -c 1 values, a (e) may be 
approximated by 



a(e) 



2d(l 



(k-rftf^r 



1 - -Tr 

a 



(k • r) d^r (r <g> r + el) 



2d(l + e) 



1 



1 



-2- + e +e 



.(A7) 
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Appendix B: Explicit results for diluted networks 

In this Appendix we present the explicit results for diluted networks. 
The nonlinear differential bulk modulus calculated using (|49[) is given by 



<P ((d + 2) 2 (Z - 2)e 2 - 2 (d 2 + 7d - 2(d + 2)2 + 4) e - 6d + 3Z) 3 

e 3 ((z - 2)e((Z - 2)e((Z - 2)e - 6) + 12) - 8)d 6 
+2e 2 (6(z - 2)(Z - 2) 2 e 4 + 3(2 - 2){Z - 2){2Z - 15)e 3 ) d 5 
+2e 2 (-3(2 - 2)(11Z - 42)e 2 + 4(142 + Z - 39)e - 36) d 5 
+e (60(z - 2)(Z - 2) 2 e 5 + 24(z - 2)(Z - 2)(5Z - 21)e 4 + 3(2 - 2)(Z(19Z - 262) + 564)e 3 ) d 4 
+e (-2((Z - 334)Z + z(137Z - 626) + 1428)e 2 + 72(5z + Z - 17)e - 216) d 4 
+2 (80(2 - 2)(Z - 2) 2 e 6 + 24(z - 2)(Z - 2)(10Z - 29)e 5 ) d 3 
+2 (12(2 - 2)(Z(19Z - 140) + 206)e 4 + ((2398 - 167Z)Z) d 3 
+2 (z(Z(68Z - 1063) + 2198) - 4636)e 3 - 9((Z - 94)Z) d 3 
+2 (z(29Z - 114) + 276)e 2 + 108(2z + Z - 7)e - 108) d 3 
+ (240(2 - 2)(Z - 2) 2 e 6 + 96(z - 2)(Z - 2)(10Z - 21)e 5 + 24(2 - 2)(Z(57Z - 274) + 282)e 4 ) d 2 
+ (9(12(51 - 5Z)Z + 2 (19(Z - 14)Z + 336) - 544)e 2 - 54(8z(Z - 2) + (Z - 28)Z + 16)e + 108(2 + 2Z)) d 2 
+ (8(5(286 - 45Z)Z + z(3Z(34Z - 231) + 742) - 1428)e 3 ) d 2 
+2 (96(2 - 2)(Z - 2) 2 e 6 + 240(z - 2)(Z(2Z - 7) + 6)e 5 ) d 
+2 (48(2 - 2)(Z - 1)(19Z - 42)e 4 + 4((1262 - 407Z)Z) d 
+2 (z(Z(204Z - 679) + 334) - 624)e 3 + 18(z(Z(19Z - 80) + 20)) d 
+2 (-2(3Z - 8)(7Z - 2))e 2 + 27(2z(Z - 8) - 7Z + 16)Ze - 27Z(2z + Z)) d 
+z (64(4e(e(e + 3) + 3) + 5)e 3 - 8Z(4e(e(4e(2e + 9) + 57) + 41) + 45)e 2 + Z 2 (4e(e + 2) + 3) 3 ) 
-8e (64e 2 (e + l) 3 - 8Ze(e(e(4e(2e + 9) + 57) + 37) + 9) + Z 2 (2e(2e(e(4e(e + 6) + 57) + 61) + 63) + 27)) 

(Bl) 



for e > e c . Below the transition curve the nonlinear differential bulk modulus vanishes. The transition curve for the 
assumption of the first order transition is given by 



18d(-2d + Z) + 3(-8d(4 + d(7 + d)) + (12 + d{29 + 7d))Z)e 
-4(16 + d(80 + d(8l + d{20 + d-2Z)- 21Z) - 48Z) - 28Z)e 2 
v +(2 + rf) 2 (-8(4 + d{7 + d)) + (28 + d(19 + d))Z)e 3 + 2(2 + d) 4 (-2 + Z)e 4 J 

( 9(-2d + Z) - 3 (4 + 2,\d + 13d 2 - 8(2 + d)Z) e - 2(2 + d)(2(8 + 5d(4 + d)) — 11(2 + d)Z)e 2 S 
-(2 + d) 2 (20 + 2hd + 3d 2 - 8(2 + d)Z) e 3 + (2 + d) A {-2 + Z)e 4 



The transition curve for the assumption of the second order transition is given by 



(d + 2) 6 (Z - 2) 2 e 6 - 6(d + 2) 4 (Z - 2) (d 2 + 7d- 2{d + 2)Z + 4) e 5 
+3(d + 2) 2 (I9(d + 2) 2 Z 2 - 2{d + 2)(d(lld + 65) + 38)Z + 4(d(d + i)(d(d + 13) + 17) + 16)) e 4 
+ ((d + 2) 2 {d(d + 163) + 244)Z 2 - 2{d + 2){d(d(d(2d + 163) + 873) + 1114) + 296)Z) e 3 
+ (4(d(d + 7) + 4)(d(d(d(d + 32) + 129) + 128) + 16)) e 3 
+9{2d - Z)(2(d(d + A)(d{d + 13) + 17) + 16) - (d + 2)(d{d + 28) + 28)Z)e 2 
+27{d{d + 7) + 4)(Z - 2d) 2 e + 27d(Z - 2d) 2 

(d + 2f{Z - 2) 2 e 6 - 6(d + 2) 4 (Z - 2) (d 2 + 7d- 2(d + 2)Z + 4) e 5 
+3(d + 2) 2 (I9(d + 2) 2 Z 2 - 2{d + 2)(d{lld + 65) + 38)Z + 4(d(d + 4){d(d + 13) + 17) + 16)) e 4 
+2(d + 2) (68(rf + 2) 2 Z 2 - (d + 2)(<f(137d + 515) + 164)Z + 2d{d(d + 5)(28d + 117) + 314) + 80) 
+9(d + 2){d(20d - 19Z + 74) - 38Z + 20)(2d - Z)e 2 + 108(d + 2)(Z - 2d) 2 e + 27(Z - 2d) 2 



(B3) 
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Appendix C: Calculation of the differential non-affinity parameter 



The only unknown term in the differential non-affinity expression in Eq. (|38|) is y ^ dv ^ e ^ V This term may be 
evaluated within the framework of the EM theory as following. Using Eqs. (|26|27I30|) one obtains 

de 



f dw n (e) 
v de 



1 



2N- 



■E 



a(e) 



(l-e- kf ) 



e(l+e) 

1 2' 

e(l+e) JJ.-jJ.jj 

,(.:) 



E(r®r+el)(l-e lk r ) ± E(r(8r+el)(l-e- ik r ) ± 

TT r T T^r 



(l-^ kr ) 



l](r®r+el)(l-e i 



j^r(r®r+el)(l-e-» k --) 
(l-e*-')E(l-e-* V ) (l-e- lkr )E(l-e- !kr ') 
|"E(rigir+eI)(l-e*r)| ^( r0r+e] [)(i_ e -<k-r; 

As before we apply the approximation of the highly coordinated lattice and get 



> . 



dv n (e) 



e(l+e) fJ-fJj 



a(e) 
e(1+<0 



,,(,-) 
Jj-fJ 



fJ+fJi. 



1 2\ 



+ 
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